function HCP_EN_predictable_behavior(FC_file, EN_dir, maxEN_iter, Nperm, test_metric, intrim_csv, outmat, ...
    restricted_csv, subj_ls, bhvr_ls, colloq_ls)

    % HCP_EN_predictable_behavior(FC_file, EN_dir, maxEN_iter, Nperm, test_metric, intrim_csv, outmat, ...
    %     restricted_csv, subj_ls, bhvr_ls, colloq_ls)
    %
    % Permute behavioral scores across subjects, considering the family structures (multi-level block 
    % permutation). Muilti-level block permutation requires the FSL PALM package.
    % Run elastic net (EN) on the permuted scores to test the predictability of the original 
    % KRR model.
    %
    % Inputs:
    %   - EN_dir
    %     The directory which contains the elastic net trained models and testing results.
    %   - maxEN_iter
    %     Maximal random seed used to split the training-test folds for performing elastic net, e.g. 400.
    %   - Nperm
    %     Number of premutations.
    %   - test_metric
    %     The accuracy metric used to perform statistical testing. Choose from 'predictive_COD' and 'corr'.
    %   - intrim_csv
    %     Full path of the intermediate csv file generated by hcp2block2 function.
    %   - outmat
    %     Full path of the output mat file storing the behaviors whose actual prediction accuracy was significantly
    %     above chance.
    %   - restricted_csv 
    %     Full path of the HCP restricted csv file, containing the family structure information. 
    %   - subj_ls
    %     Full subject list (absolute path).
    %   - bhvr_ls 
    %     List of behavioral measures for which matched AA and matched WA can be found (absolute path).
    %   - colloq_ls
    %     List of colloquial names of behavioral variables. The colloquial names should correspond to the
    %     behavioral names in "bhvr_ls".

    addpath(fullfile(getenv('CBIG_CODE_DIR'), 'external_packages', 'matlab', 'non_default_packages', 'palm', 'palm-alpha109'))
    addpath(fullfile(getenv('CBIG_CODE_DIR'), 'external_packages', 'matlab', 'non_default_packages', 'glmnet', 'glmnet_matlab'))

    iter = 40;
    
    %% load shared files
    [bhvr_nm, nbhvr] = CBIG_text2cell(bhvr_ls);
    [colloq_nm, ncolloq] = CBIG_text2cell(colloq_ls);
    if(ncolloq ~= nbhvr)
        error('Number of behavioral names is not equal to number of colloquial names.')
    end
    [subjects, nsub] = CBIG_text2cell(subj_ls);
    alpha_FDR = 0.05;
    
    %% create permutations
    % 1. use `hcp2blocks.m` to generate block definitions
    % 2. use `palm_tree.m` to create exchangeable tree
    % 3. use `palm_permtree.m` to generate permutation indices
    idx_perm_out = fullfile(EN_dir, 'Pset.mat');
    if(~exist(idx_perm_out, 'file'))
        B = hcp2blocks(restricted_csv, intrim_csv, false, dlmread(subj_ls));
        Ptree = palm_tree(B);
        Pset = palm_permtree(Ptree, Nperm+1, false, true);
        save(idx_perm_out, 'Pset', 'B')
    else
        load(idx_perm_out)
    end
    
    %% load and reshape RSFC 
    load(FC_file)
    features = reshape3d_to_2d(corr_mat);

    %% Read covariates which need to be regressed out from RSFC
    cov_X = load(fullfile(EN_dir, 'cov_X_58behaviors.mat'));

    %% calculate permuted elastic net accuracies
    metrics = {'corr','COD','predictive_COD','MAE','MAE_norm','MSE','MSE_norm'};
    parfor b = 1:nbhvr
        EN_perm(features, cov_X.covariates, b, EN_dir, maxEN_iter, Nperm, Pset, bhvr_nm, metrics)
    end

    %% calculate p value for each behavior
    p_perm = zeros(nbhvr, 1);
    avg_stats = [];
    avg_null_stats = []
    for b = 1:nbhvr
        curr_avg_stats = [];
        curr_avg_null_stats = [];
        for i = 1:maxKRR_iter
            opt_fname = fullfile(EN_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'optimal_acc', ...
                [bhvr_nm{b} '_final_acc.mat']);
            if(~exist(opt_fname, 'file'))
                continue
            end
            opt = load(opt_fname);
            Nfolds = length(opt.original_statistics);
            orig_stats = zeros(Nfolds, 1);
            for f = 1:Nfolds
                orig_stats(f) = opt.optimal_statistics{f}.(test_metric);
            end

            acc_out = fullfile(EN_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'perm.mat');
            load(acc_out)
    
            curr_avg_stats = cat(1, curr_avg_stats, mean(orig_stats, 1));
            curr_avg_null_stats = cat(1, curr_avg_null_stats, mean(stats_perm.(test_metric), 1));
        end
    
        avg_stats = cat(2, avg_stats, curr_avg_stats);
        avg_null_stats = cat(3, avg_null_stats, curr_avg_null_stats);
        p_perm(b) = length(find( mean(curr_avg_null_stats,1) > mean(curr_avg_stats,1) )) ./ Nperm;
    end

    %% multiple comparisons correction
    H = FDR(p_perm, alpha_FDR);
    sig_perm_idx = sort(H);
    sig_behaviors = bhvr_nm(sig_perm_idx);
    
    switch test_metric
    case 'predictive_COD'
        sig_COD = sig_perm_idx;
        p_COD_perm = p_perm;
        save(outmat, 'p_COD_perm', 'H', 'sig_COD', 'sig_behaviors', 'avg_stats', 'avg_null_stats');
    case 'corr'
        sig_corr = sig_perm_idx;
        p_corr_perm = p_perm;
        save(outmat, 'p_corr_perm', 'H', 'sig_corr', 'sig_behaviors', 'avg_stats', 'avg_null_stats');
    end

    rmpath(fullfile(getenv('CBIG_CODE_DIR'), 'external_packages', 'matlab', 'non_default_packages', 'palm', 'palm-alpha109'))
    rmpath(fullfile(getenv('CBIG_CODE_DIR'), 'external_packages', 'matlab', 'non_default_packages', 'glmnet', 'glmnet_matlab'))
    
end


function EN_perm(features, cov_X, b, EN_dir, maxEN_iter, Nperm, Pset, bhvr_nm, metrics)

    % EN_perm(features, cov_X, b, EN_dir, maxEN_iter, Nperm, Pset, bhvr_nm, metrics)
    %
    % 

    for i = 1:maxEN_iter
        load(fullfile(EN_dir, ['randseed_' num2str(i)], bhvr_nm{b}, ...
            ['no_relative_10_fold_sub_list_' bhvr_nm{b} '.mat']));
        Nfolds = length(sub_fold);

        acc_out = fullfile(EN_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'perm.mat');
        if(exist(acc_out, 'file'))
            continue
        end

        for f = 1:Nfolds
            test_ind = sub_fold(f).fold_index==1;
            train_ind = ~test_ind;

            %% load cross-validatedly regressed behavioral scores
            y_reg = load(fullfile(EN_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'y', ...
                ['fold_' num2str(f)], ['y_regress_' bhvr_nm{b} '.mat']));

            %% split training vs test data for RSFC
            feat_train = features(:, train_ind)';
            feat_test = features(:, test_ind)';

            %% do confound regression from features, if necessary
            if(~isempty(cov_X) && ~strcmpi(cov_X, 'none'))
                [feat_train, beta] = CBIG_regress_X_from_y_train(feat_train, ...
                    cov_X(train_ind, :));
                beta_pre = load(fullfile(EN_dir, ['randseed_' num2str(i)], bhvr_nm{b}, ...
                    'params', ['fold_' num2str(f)], 'feature_regress_beta.mat'));
                if(max(abs(beta(:) - beta_pre.beta(:))) > 1e-8)
                    error('[Regression from RSFC]: beta differred from original elastic net results')
                end

                feat_test = CBIG_regress_X_from_y_test(feat_test, ...
                    cov_X(test_ind, :), beta);
            end

            %% load optimal parameters selected in inner-loops
            params_file = fullfile(EN_dir, ['randseed_' num2str(i)], bhvr_nm{b}, 'params', ...
                ['fold_' num2str(f)], ['selected_parameters_' bhvr_nm{b} '.mat']);
            opt_params = load(params_file)

            %% multilevel block permute y, run elastic net based on permuted y
            for p = 2:Nperm+1
                rng('default')
                rng(1)
                opts.lambda = opt_params.curr_lambda; % should take best lambda value from innerloop cv here
                opts.alpha = opt_params.curr_alpha;

                % permute y
                y_perm = y_reg.y_resid(Pset(:,p));
                y_train = y_perm(train_ind);
                y_test = y_perm(test_ind);

                fit=glmnet(feat_train, y_train, [], opts);
                y_p = glmnetPredict(fit, feat_test);

                for k = 1:length(metrics)
                    stats_perm.(metrics{k})(f,p-1) = ...
                        CBIG_compute_prediction_acc_and_loss(y_p, y_test, metrics{k}, y_train);
                end
            end
        end
        save(acc_out, 'stats_perm')
    end
    
end

function out = reshape3d_to_2d(features)
    % reshapes a #ROI x #ROI x #subjects matrix into
    % #ROI x #subjects by extracting the lower triangle
    % of the correlation
    temp = ones(size(features,1), size(features,2));
    tril_ind = tril(temp, -1);
    features_reshaped = reshape(features, size(features,1)*size(features,2), size(features, 3));
    out = features_reshaped(tril_ind==1, :);
end